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ABSTRACT 


Pulse Detonation Engines (PDE) are seen to be the next generation propulsion 
systems due to enhanced thermodynamic efficiencies based on the Humphrey cycle. One 
of the limitations in fielding practical designs has been attributed to tube diameters not 
exceeding 5 inches as the shock wave takes a long run distance for transition to 
detonation, thus potentially affecting specific thrust. Novel methods via imploding 
detonations were investigated to remove such limitations. During the study, a practical 
computational cell size was first determined so as to capture the required physics for 
transient detonation wave propagation using a Hydrogen-Air test case. Through a grid 
sensitivity analysis, one-quarter of the induction length was found sufficient to capture 
the experimentally observed initial wave transients. Test case models utilizing axi- 
symmetric head-on implosions were studied in order to understand how the implosion 
process reinforces a detonation wave as it expands. This in effect creates localized 
overdriven regions, which maintains the transition process to full detonation. A 
parametric study was also performed to determine the extent of diameter increase and it 
was found that the detonations could be supported with no change in run distance even 
when the tube diameter exceeds 5 inches. 
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I. INTRODUCTION 


A. MOTIVATION 

Pulse Detonation Engines (PDEs) have tremendous potential to replace current 
state-of-the art propulsion systems, which can operate up to supersonic cruise speeds for 
missiles, fighter jets, unmanned aerial vehicles, etc. Much of this potential stems from the 
fact that PDEs operates on the detonation mode of combustion, which has a higher 
thermodynamic efficiency than current propulsion systems based upon constant pressure. 
However, no practical design has been fielded to date and one key limitation is due to the 
fact that detonation can only be supported in main detonation chambers of up to five 
inches in diameter without requiring unreasonably long tube lengths. This impairs the 
specific thrust substantially as multiple tubes are required to generate thrust required by a 
typical cruise missile. As a result, thrust density decreases and potential cross flow 
mixing problems from the efflux of five separate tubes reduce the efficiency even further. 
To increase the diameter of the main tube and yet maintain the tube length, simulation 
studies involving novel means using shock-focusing techniques via toroidal wave 
implosion first investigated by Jackson et al. [1] is conducted. This reinforces the wave 
propagation process and effectively shortens the deflagration-to-detonation run distance 
such that a detonation wave can be supported for larger diameter tubes. 

B. BACKGROUND 

Pulse Detonation Engines (PDEs) operates on the detonation mode of 
combustion, which can be best approximated by the Humphrey (Constant-Volume) cycle 
in contrast to the Brayton (Constant-Pressure) cycle for gas turbine engines, as illustrated 
in Figure 1. The amount of work done (Pressure-Volume) by the Humphrey cycle as heat 
energy is introduced under constant volume conditions from 2 to 3 is much greater as 
compared to the Brayton cycle under constant pressure conditions, resulting in higher 
thermodynamic efficiencies. 
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Figure 1. Pressure-volume plot based on Humphrey cycle [From 2] 

In order to create detonation shock waves that can propagate at supersonic speeds, 
a finite run length distance is required for such transitions from deflagration-to- 
detonation known as a DDT process for a typical confined tube configuration. As tube 
diameters become larger, a longer DDT run distance is required in order to support the 
transition to a fully developed detonation wave. This limitation results in multiple tube 
configurations for a practical cruise missile design. 

A novel method is being proposed using shock focusing techniques via toroidal 
wave implosion such that the shock wave propagation in the main shock tube is 
reinforced and consequently, the resulting pressures generated are much higher than the 
Chapman-Jouguet (C-J) conditions defined for a fully-developed steady state wave. The 
extreme high pressure enviro nm ent favors much higher reaction rates for the fuel- 
oxidizer reactant mixtures and this will allow for a near direct generation of a detonation 
(without the need for DDT) as the fiiel-air mixtures react violently in order to support the 
shock front. 

In order to explore such implosion concepts accurately, it is imperative to capture 
the full physics and chemical kinetics behind such complex and highly transient fluid 
flows governing the shock wave behavior followed by chemical reactions based on a 
reduced set of chemical reaction mechanisms. Simulations of such complex flow field 
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behavior are only achievable via a practical computational cell mesh topology that is 
sufficiently small enough on the order of the induction length to describe the reaction 
zone behind the propagating shock front. Validation of such simulation studies based on 
well known reduced kinetic models for Hydrogen-Air mixtures can be conducted via 
comparison with experimental data available in the literature. 

The required computational cell size for such high fidelity simulations can then 
be further extended to generic hydrocarbon fuel-air mixtures such as acetylene-air, 
ethylene-air, propane-air and even JP 10-air mixtures based on the required length of the 
induction zone supporting the shock front. Utilizing the practical computational cell size, 
further analysis for a simulation case model based on an axi-symmetric head-on toroidal 
implosion can be conducted and analyzed for similar configuration types. 

Parametric studies are then performed to investigate the extent of the increase in 
tube diameters for a fixed practical propulsion tube length based on the test case model. 

C. OBJECTIVES 

The research objective is to enhance current PDE applications by carrying out 
simulation studies to investigate the extent of possible combustor diameter increase using 
shock focusing techniques via toroidal implosions with modified PDE tube 
configurations without introducing the need for longer tube lengths. 

Simulation studies were conducted on the CFD-I-+ computational platform and 
carried out via the following order: 

i) Determination of a practical computational cell size for high fidelity 
capturing of detonation wave propagation in a simple PDE tube 
configuration based on well known Hydrogen-Air reduced reaction 
models 

ii) Simulation and analysis of head-on axi-symmetric toroidal implosions for 
Hydrogen-Air mixture based on the required computational cell size 
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iii) Performance of parametric studies used to evaluate the extent of the 
diameter increase for a fixed practical tube length using toroidal implosion 
based on an enhanced PDE design configuration 
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II. THEORY 


A. PULSE DETONATION ENGINES 

Pulse Detonation Engines are essentially intermittent propulsion systems where 
thrust is provided by repetitive triggering of detonation waves propagating at supersonic 
speeds in a long tube enclosed at the head end. The process takes place under nearly 
constant volume conditions and its operating cycle can be effectively described by three 
distinct phases. The first phase consists of introducing a reactive component such as 
hydrogen fuel and an oxidizing component such as air into the main combustion 
chamber. Initiation of the fuel-air mixture near the enclosed portion of the chamber 
occurs in the second phase via an energy sparking source. A high pressure detonation 
wave as a result of the initiation ensues in the third phase. The remaining fuel-air mixture 
in the tube is compressed in a thin reaction zone by the traversing detonation wave and 
reacts violently under extreme high pressure conditions to feed the leading edge of the 
shock front followed by an expansion zone of detonation products. The last phase 
consists of effectively purging the gaseous products out of the main chamber and the 
entire cycle repeats. Thrust is produced as the gaseous products exit out of the main 
propulsion tube. Figures 2 to 4 describes each phase of a PDE cycle. 



Figure 2. PDE cycle - fill phase 
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Figure 3. PDF cycle - fire phase 
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Figure 4. PDE cycle - purge phase 


B. DETONATION THEORY 

Detonation or deflagration by itself is an exothermic reaction caused by transition 
of unstable reactive molecular reactions to a stable state. The amount of energy required 
to trigger such reactions is known as the Activation Energy. When such energy levels are 
reached, the original reactant molecules come apart and re-arrange themselves back into a 
more stable state producing shock waves as a result of a sudden sharp increase in 
thermodynamic properties due to the exothermic reaction. 

Combustion of fuel-air mixtures produces either a deflagration or detonation 
wave. Properties such as pressure (P), temperature (T) and density (p) of the reactant and 
products across the combustion wave front change discontinuously. Figure 5 illustrates 
the change across such properties for a constant area duct where the flame front is 
envisioned as stationary. 
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Figure 5. One-dimensional combustion wave model [From 3] 


In order to relate properties across a one-dimensional combustion wave front, the 
corresponding mass conservation, momentum conservation, energy conservation and 
equation of state must be invoked. For the state equation, the ideal gas law serves as a 
good approximation to model the properties of the gases during combustion. These 
relations assume that the combustion wave is steady state and the process across the wave 
front is adiabatic and remain in chemical and thermodynamic equilibrium [3]. 

Mass Conservation: 

p^u^=p^u^ ( 1 ) 

Momentum Conservation: 

+ /?[Mj = Pj + PjUj (2) 

Energy Conservation: 

= (3) 

Equation of State (Ideal Gas Law): 

P = pRT (4) 

Specific heat relation: 

C,=—R ( 5 ) 

Y-\ 


1 



where: 


R = Gas Constant 

q = Specific heat energy added to the system via the combustion process 
Cp = Specific heat at constant pressure 
y = Ratio of specific heats 

In order to gain further insight into all possible states that the burned gas products 
can exist, a Hugoniot relation must be determined. A Hugoniot curve is basically a plot of 

pressure (P 2 ) versus specific volume (—) corresponding to the burned gas products for 

Pi 

a fixed heat release mass, q and the relation can be derived by combining equations (2), 
(3) and (4) to obtain the given Hugoniot relation. 

Hugoniot Relation: 


r {Pi ^ 1 ' 
r-APi Ply 




( 6 ) 


The slope of the Hugoniot relation indicates the velocity of the detonation wave. The 
Rayleigh Line relation is obtained by combining equations (1) and (2). 

Rayleigh Line Relation: 


2 2 

A Ml 



Pi Pi 


( 7 ) 


The Hugoniot curve is illustrated in Figure 6. From the curve, the origin A 

represents the pressure (Pj) and specific volume (—) of the unbumed reactants. The 

Pi 

two tangent points intersecting origin A are defined as upper and lower Chapman-Jouget 
(C-J) points. From the Hugoniot curve, the C-J points can be obtained via the slope of the 
Hugoniot curve as follows: 
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C-J relation: 


dP. 


iP2-P,) 


(V ] 

f—--1 

[/plj 

K.P 2 Ay 


( 8 ) 


The C-J relation ean also be obtained via differentiation of equation (6) with respect to 
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( 9 ) 


The speed of sound, ci for the burned products is given by the following relation: 

C2=4^^2 ( 10 ) 

Equating (8) and (9) along with the Rayleigh-Line relation in (7) gives the following: 


yPi 


Pi 


Cl 


( 11 ) 


This implies that the particle velocity of the gases leaving the detonation wave are 
traveling at sonic velocity. In a shock-frame of reference, this means that the product 
gases leaving the detonation wave are leaving at the sonic velocity of the products. 
Therefore equilibrating the amount of energy that reaches the shock front from the 
chemical reactions and thus allowing a sustained steady-state detonation wave to 
propagate at what is known as the Chapman-Jouguet velocity. Should this not be the 
case, detonation would undergo an unsteady process that will most likely asymptote to 
the C-J condition. 

Five distinct regions can be inhered from the Hugioniot curve, however, not all 
regions are physically valid. Table 1 describes the validity of the states in these regions. 


9 




Figure 6 . Hugoniot curve for fuel-air mixtures [From 3] 


Table 1. Physical validity of states on the Hugoniot curve 


Region 

Remarks 

I 

Region I is defined as transient states where the combustion waves propagate 

as strong overdriven detonations above the upper C-J point. Under such 

conditions, the lead shock front travels at a much faster rate than the product 

gases resulting in a larger induction zone. The induction zone will eventually 

increase to such an extent such that any heat release generated by the 

combustion of the fuel-oxidizer reactants will not affect the lead shock front, 

thus slowing it down. As such, strong detonation waves will always decay 

back to the upper C-J point. 

II 

Region II is defined as transient states where the combustion waves propagate 
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Region 

Remarks 


as weak detonations below the upper C-J point. Under such conditions, the 

velocity of the gas products travels at a faster rate than the lead shock front and 

as such the spatial extent of the induction zone decreases. The same heat 

release energy generated from the combustion reaction starts to have an 

increasing effect on the lead shock front due to the decreasing induction zone 

length and the additional heat release supplied starts to accelerate the lead 

shock front away from the reaction zone back towards the upper C-J point. 

III 

Region III represents states where the combustion waves propagate as weak 

deflagrations. They exist as expansion or rarefaction waves where the density 

of the burned products is lower than that corresponding to the reactants. They 

are observed experimentally during the initial deflagration-to-detonation 

process at ignition and the gas velocity relative to the wave front is accelerated 

from lower to higher sub-sonic speeds. 

IV 

Region IV represents states where the combustion waves propagate as strong 

deflagration waves. These are physically impossible states where the gas 

velocity relative to the wave front must be accelerated from sub-sonic to super¬ 
sonic speeds, which are forbidden for a wave structure propagating in a 

constant area duct. 

V 

The solution states for region V is mathematically impossible. The solution to 

the Rayleigh-Line equation (7) results in an imaginary number for ui. 


From the analysis given in Table 1, the physically validity of regions 1 and 3 are 
most likely observed during a detonation process and this has been confirmed in shock 
tube experiments as well [3]. 

C. DEFLAGRATION-TO-DETONATION TRANSITION 

A detonation wave requires a finite run distance for the shock wave to transit from 


a deflagration to a detonation wave known as DDT. During initiation of a fuel-oxidizer 
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mixture in a PDE tube, a deflagration leading shock wave ensues and propagates at a 
wave speed characteristic of the reactant temperature conditions followed by a series of 
compression waves produced by the combustion process. The compression waves further 
heats up the region behind the leading shock wave creating localized high temperature 
regions and consequently increase the velocity of the compression waves. The 
compression waves eventually catch up with the leading wave front coalescing into a 
detonation wave [3]. The tube confinement as well as the turbulence generated by the 
motion of the product gases from the detonation shock front leads to the onset of “an 
explosion in an explosion,” which produces an overdriven forward shock wave into the 
unbumed reactants known as superdetonation and an opposite shock wave moving into 
the gas products known as retonation with transverse oscillating shocks in between. This 
finally leads to a self-sustained steady state detonation wave propagating at C-J velocity. 
Figure 7 illustrates this phenomenon. 



Figure 7. Transition process to superdetonation [From 3] 


The DDT run distance is impractically long for current applications and therefore 
obstacles such as swept rams or spirals are used in PDE configurations to shorten the 
DDT length via turbulent flow regimes to promote higher chemical reaction kinetics. 
However, generation of high pressure losses as a consequence of such obstacles degrades 
the propulsive efficiency tremendously. 


12 



D. DETONATION WAVE STRUCTURE 


Sustainment of the detonation wave can be described simply by the one¬ 
dimensional model developed independently by Zeldovich, von Neumann and Doring 
[3]. This is known as the ZND model illustrated in Figure 8. The model consists of a 
shock wave moving at detonation velocity followed by a rarefaction region. The shock 
wave region is extremely thin in the order of a few mean free paths thick of molecules; 
therefore limited reactions occur in this region. Across the shock region, there is a sharp 
almost discontinuous increase of pressure, temperature and density. The rarefaction 
region is comparatively much thicker usually on the order of 1 cm and consists of an 
induction and reaction zone. The thermodynamic properties in the induction zone 
immediately behind the shock front remain almost constant. Reaction rates as a function 
of temperature begin to increase slowly. Once the induction period is completed, 
temperatures start to climb sharply with extremely high heat release energies caused by 
the violent reactions of the fuel-oxidizer mixture defined by the reaction zone, which 
sustains the propagation of the shock wave. The gasdynamic properties reach their 
equilibrium states at the end of the reaction zone. Therefore, the induction period or 
induction zone is extremely important in order to capture the correct thermodynamic 
properties of temperature, pressure and density for an accurate simulation of a detonation 
wave. 
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Figure 8. Thermodynamic properties of the ZND model [From 3] 


E. DETONATION CELL WIDTH 

Detonation cell widths are multi-dimensional detonation wave effects 
experimentally measured by soot foil imprint techniques where a detonation wave created 
by a typical fuel-air mixture leaves a fish scale pattern signature on soot coated aluminum 
sheet lined on the inner surface of a shock tube. Figure 9 illustrates a typical soot foil 
imprint produced by a detonating shock wave from JPlO-Air mixture where JPIO is a 
liquid hydrocarbon fuel (CioHie). The non-planar detonation shock front induced by the 
energy release arising from the combustion of the fuel-oxidizer mixture is a resultant of 
the Mach stem, incident shock and reflected shock, which is a necessary condition for a 
self-sustained detonation front. The interaction of the three waves produces a shear 
discontinuity known as the triple point. Figure 10 illustrates such a detonation front 
structure. As a result, the waves propagate down the tube leaving the characteristic fish 
scale pattern signature. The transverse distance between each fish scale is known as the 
cell width, which is regarded as a fundamental parameter to characterize the detonability 
of a gaseous fuel-air mixture. A more reactive fuel such as hydrogen as compared to a 
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typical hydrocarbon fuel such as ethylene (C 2 H 4 ) will correspondingly have a smaller 
detonation cell width. Similarly, reaction of hydrogen in the presence of a less reactive 
oxidizer such as air with respect to pure oxygen will have larger cell widths. 



Figure 9 Soot foil imprint produced by JPlO-Air detonation wave [From 4] 
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Figure 10. Detonation front structure with triple point [From 4] 
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F. 


REDUCED CHEMICAL KINETIC REACTION MODELS 


Chemical kinetics governs the rate at which the species concentration of a fuel- 
oxidizer mixture is transformed into its product species and vice versa depending on 
pressure and temperature reaction conditions as well as complexity of the involved 
species. The dependence of the reaction rate constant, k as a function of temperature, T 
can be inferred from the Arrhenius Equation [3]. 

Arrhenius Equation: 


f 

k = AT" exp 

V 



RT ) 


( 12 ) 


The pre-exponential parameter A describes the frequency of collisions for the species that 
result in a reaction as well as the preferred orientation of the molecular species that favor 
a reaction. The exponent n expresses the temperature dependence of the pre-exponential 
parameter. Activation energy, Ea defines the minimum amount of energy required for a 
reaction to occur and R is the universal gas constant. 

While it is desired to include a detailed chemical reaction model governing a 
reactive fuel-oxidizer mixture, which can involve up to hundreds of reactions and species 
under all operating conditions, it is computationally demanding especially for complex 
flows involving cellular detonation wave structures even for super-computers of today. 
To ensure that computational times required for CFD simulations to be tractable under 
reasonable time frames, it is imperative to reduce detailed chemical mechanisms into 
simpler models where only certain chemical reactions are significant for example during 
a detonation process and yet able to capture the physics of the shock wave propagation 
with correct heat release rates in the reaction zone. This can be done via approximation of 
individual reactions to be in equilibrium or some of the species concentration can be 
assumed to be in quasi-steady state [5]. 


G. EQUIVALENCE RATIO OF FUEL-OXIDIZER MIXTURE 

Typically, combustion of fuel-oxidizer mixtures is specified in terms of the 
equivalence ratio, c|). The equivalence ratio is independent of the mass, m or the number 
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of moles, n that is used to eharacterize the proportion of the fuel to oxidizer mixture. It is 
defined as the ratio of the fuel-oxidizer ratio to the stoiehiometrie fuel-oxidizer ratio and 
can be mathematically represented as follows: 


<!> = 


^fuel I ^oxidizer 


^ fuel j ^oxidizer 


\ fuel / oxidizer /stoichiometric 


fuel j ^oxidizer 


stoichiometric 


(13) 


The composition and thermodynamic properties of fuel-oxidizer mixtures during 
combustion/detonation are dependent on the equivalence ratio and consequently the 
detonation cell width if the wave propagates at supersonic speeds. A ratio greater than 1 
would imply excess fuel or fuel rich mixtures than would be required for a stoichiometric 
reaction and conversely a ratio less than 1 would imply excess oxidizer or fuel lean 
mixtures. 


H. COMPUTATIONAL NUMERICAL SCHEMES 

In computational fluid dynamics, there are primarily two main schemes, explicit 
and implicit schemes implemented in order to discretize the resulting differential 
equations and resolve the required flow field. Explicit or forward differencing schemes 
suffer from conditional stability issues where a solution may blow up if the specified time 
step is too large and/or the computational cell size is too small. Implicit schemes are 
inherently stable as the forward solution is coupled together with the previous solution 
and hence will not generate an output that will increase to infinitely large values. 
Schemes of such types will be employed as the main scheme in study for simulation of 
detonation waves in PDE configurations. Figure 11 show a computational molecule based 
on an implicit scheme where p defines the spatial grid points in space and q defines the 
temporal grid points on an x-t diagram. 
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Figure 11. Computational molecule based on implicit scheme 


For a mathematical representation of an implicit scheme, consider a simple one¬ 
dimensional scalar equation: 


du du 

— + a — = 0 
dx dt 


(14) 


Discretization of the above equation based on an implicit scheme results in the following: 


a+\ 


-ul+- 


aAt ( 
\u 


Ax 


q+\ 

p+l 


-< l )=0 


( 15 ) 


-is also known as the Courant-Fredrichs-Lewy or CFL number. Numerical solution of 

Ax 

the resulting flow field can be specified either in terms of a fixed time step or the CFL 
number where the time step is computed based on the computational mesh size. For a 
sufficiently small fixed time step, specification of the CFL number will result in faster 
numerical computations [6]. 
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III. SIMULATION SETUP (GRID SENSITIVITY) 


A. WAVE IMPLOSION 

Wave implosion is a technique that focuses a shock front to generate a region of 
high pressure and temperature gases due to adiabatic compression of gas products as it 
flows into a region of decreasing area [1]. As a result, extremely high energy density 
localized regions are formed due to interacting compression waves and this raises the 
pressure of a typical fuel-oxidizer mixture above the C-J pressure, resulting in an 
overdriven detonation wave. Figure 12 illustrates an example of a collapsing toroidal 
wave front implosion. This concept is being explored to increase specific thrust of PDFs 
by harnessing the extremely high pressures as a result of the overdriven imploding shock 
front to increase the maximum diameter for a DDT process assuming a reasonable run 
distance required for transition to detonation. A commercial Computational Fluid 
Dynamics (CFD) platform used to perform simulation studies on the implosion concepts 
for PDFs is CFD-i-i- from Metacomp Technologies. It is widely used by commercial 
companies such as General Flectric and Air Force Research Laboratories (AFRL) as a 
computational tool for development on proof of concept PDFs and also supported by 
Office of Naval Research (ONR). 
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B. HIGH FIDELITY SIMULATIONS (GRID SENSITIVITY ANALYSIS) 


In order to explore implosion concepts, high fidelity computational fluid dynamic 
(CFD) models based on reduced chemical kinetic mechanisms are pertinent for accurate 
and robust simulation to describe and analyze the physics of detonation wave propagation 
in pulse detonation engines. High fidelity models are sensitive to the induction length for 
a specific fuel-oxidizer mixture such that the transient flow kinetics and shock physics 
can be accurately captured as shown by the ZND model in the previous section. 

Using a test case model for Hydrogen-Air (H 2 -Air) mixtures, a grid sensitivity 
analysis was performed on a typical PDE tube configuration using a kinetic reaction 
model based on a reduced set of 18 reactions and 9 species [7] and [ 8 ] highlighted in 
Table 2. 


Table 2. Reduced chemical kinetic mechanism for H 2 -Air mixture 


Number 

Reaction 

A** 

n'* 

F ^ 

1 

H2 + O 2 ^ 20H 

0.17ell 

0 

0 . 20 15e9 

2 

O 2 + H ^ 20H 

0.142el2 

0 

0.6866e8 

3 

H 2 + OH ^ H 2 O + H 

0.316e5 

1.8 

0.12686e8 

4 

H 2 + 0 ^ OH + H 

0.207el2 

0 

0.57568e8 

5 

20H ^ H 2 O + 0 

0.55ell 

0 

0.29307e8 

6 

OH + H + M ^ H 2 O + M 

0.2211el7 

-2 

0 

7 

H + H + M ^ H 2 

0.653el2 

-1 

0 

8 

O 2 + H + M ^ HO 2 + M 

0.32el3 

-1 

0 

9 

OH + HO 2 ^02 + H 2 O 

0.5ell 

0 

0.4186e7 

10 

H + HO 2 ^ H 2 + O 2 

0.253ell 

0 

0.29307e7 

11 

H + HO 2 ^ 20H 

0.199el2 

0 

0.7536e7 

12 

0 + HO 2 ^ O 2 + OH 

0.5ell 

0 

0.4186e7 

13 

2 HO 2 ^ O 2 + H 2 O 2 

0.199el0 

0 

0 

14 

H 2 + HO 2 ^ H + H 2 O 2 

0.301e9 

0 

0.7829e8 

15 

OH + H 2 O 2 ^ H + H 2 O 2 

0.102ell 

0 

0.7954e7 

16 

H + H 2 O 2 ^ H 2 O + OH 

0.5el2 

0 

0.4186e8 

17 

0 + H 2 O 2 ^ OH + HO 2 

0.199ell 

0 

0.247e8 

18 

H 2 O 2 + M ^ 20H + M 

0.121el5 

0 

0.19049e9 


in units of kmol,m^, s, K and J/kmol 
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Nitrogen and carbon chemistry reactions were excluded in the reaction model so 
as to optimize the required computational effort. Only slight differences in CJ detonation 
wave speeds and temperatures for typical hydrocarbon fuel-air mixtures on the order of 
6% to 9% respectively [9] were observed as compared to the full kinetic model 
comprising the full nitrogen and carbon chemistry for Jet A/Air mixture, as shown in 
Figures 13 and 14, which was deemed as acceptable. 
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Figure 13. Comparison of detonation velocities for Jet A/Air mixture [From 9] 
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Figure 14. Comparison of C-J temperatures for Jet A/Air mixtures [From 9] 
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C. DETONATION CELL WIDTH OE H 2 -AIR FOR DIFFERENT 
EQUIVALENCE RATIOS 

Determination of the computational mesh size for the model was focused mainly 
on stoichiometric H 2 -Air mixtures. This is because it is well defined in literature that 
experimental results at which detonation cell width occurs at a minimum corresponds to 
= 1 as indicated in Figure 15 referenced from [10], [11], [12] and [13] for a H 2 -Air 
mixture. The induction length is in turn proportional to the detonation cell width and the 
relation is typically a linear one with a proportionality constant for different fuel-oxidizer 
mixtures. Therefore, this would present the baseline mesh size required for accurate 
simulation of detonation shock physics. 
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Figure 15. Graph of detonation cell width vs equivalence ratio for H 2 -Air 

mixture [From 14] 
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D. 


PDE CONFIGURATION 


For the grid sensitivity analysis, a 2-D axisymmetric cylinder model was 
generated on the Solid Works platform with the dimensions outlined in Figure 16. A 
detonation wave consisting of the detonation products was assumed to be initiated at CJ 
Pressure and CJ Temperature conditions characteristic of the fuel-oxidizer mixture as part 
of the initial conditions. Mass fraction of the detonation products are obtained from 
NASA Chemical Equilibrium Analysis (CEA) [15] and included in Appendix A. The 
reactant initial conditions were assumed to be at 1 atmosphere and 298 K, which would 
mimic experimental conditions for further validation. 



Figure 16. 2-dimensional axisymmetric cylinder model 


The model was gridded on the Multipurpose Intelligent Meshing Environment 
(MIME) platform for an unstructured mesh (triangles) topology, which was optimized as 
a pre-processor for CFD++. It would have been ideal to resolve the flow field in terms of 
a structured mesh topology (quadrilaterals) so as to minimize the computational 
workload; however MIME did not have the capability to generate structured grids. An 
illustration of the cylinder model for unstructured mesh topology can be found in Figure 
17. 
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Figure 17. 2-D axisymmetric model (Unstructured mesh size) 


For the implosion study, a modified PDE configuration was conceived based on a 
2-D cylinder model as depicted in Figure 18 with an axi-symmetric axis. This model 
utilizes the concept of implosion with an assumption of a starting detonation wave 
propagating at the inlet. The lead shock front of the detonation wave should be sustained 
by the implosion concept and will not decay to a rarefaction wave as it expands into the 
tube space. Further parametric studies are then conducted by increasing the radius of the 
model until the detonation wave fails, thereby compromising the DDT length that is 
required for such geometric configurations. 
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E. SIMULATION PARAMETERS 

The simulation parameters used to perform the grid sensitivity analysis is 
summarized in Table 3. The unsteady compressible Euler equation (no viscosity terms) 
was chosen as the main equation set governing the numerical computations. The 
parameter of interest was to evaluate the thermodynamic parameters and wave speed 
along the symmetry axis so as to capture the detonation wave physics rather than near the 
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walls (where viscous effects become significant). A Total Variation Diminishing (TVD) 
limiter scheme based on the minmod function was used in order to prevent spurious 
oscillation of the resolved flow field in the vicinity of strong gradients, which represent 
shocks. Near flow field regions where shocks are encountered, they switch the spatial 
discretization scheme to a first order accurate method while second-order accurate 
numerical schemes are computed for the rest of the flow field. 


Table 3. Simulation setup on CFD++ platform 


Equation Set Type 

i) 

Unsteady Compressible Euler 

ii) 

Equation of State: Ideal Gas 

Riemann Solver 

i) 

Minimum Dissipation: LHS & RHS 

ii) 

Activate Pressure Switch: Supersonic 

hi) 

Activate Pressure Gradient Normal 

Boundary Conditions (See Figure 17) 

i) 

Top: Multi-species inviscid surface tangency (Wall) 

ii) 

Left: Multi-species inviscid surface tangency (Wall) 

hi) 

Right: Multi-species inviscid surface tangency (Wall) 

iv) 

Bottom: X-axis symmetry 

Boundary Conditions (See Figure 18) 

i) 

Inlet: Specified by inlet stagnation temperature, pressure and 
detonation products 

ii) 

Wall: Multi-species inviscid surface tangency 

hi) 

Symmetry: X-axis symmetry 

Time Integration 

i) 

Implicit 

Spatial Discretization 

i) 

2^^ Order Accuracy in Space 

ii) 

Dimension of polynomial: 2-D axisymmetric 

hi) 

Axis stability enhancement: Yes 

iv) 

Types of Total Variation Diminishing (TVD) Limiter: Minmod 
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IV. DISCUSSION OF RESULTS AND ANALYSIS 


A. INITIAL COMPARISON OF DIFFERENT GRID SIZES 

A grid size from 0.25 mm to 1.0 mm was first considered in the grid sensitivity 
analysis so as to determine a computational mesh size that is sufficiently small enough to 
capture the physics of the detonation wave. The ZND induction length for detonation of 
stoichometric H 2 -Air mixtures (0.23 mm) was used as a baseline comparison [16]. 

As outlined previously in the PDE configuration description, a detonation wave 
was assumed to propagate at the start with conditions pertaining to C-J pressure and 
temperature obtained from NASA CEA code. This is outlined in Table 4. (Note that the 
initial reactants were assumed to be at Pressure = 1 atm and Temperature = 298K). This 
would mimic experimental conditions for purpose of validation. The single step reaction 
equation is indicated in Table 5. 


Table 4. C-J parameters for stoichometric Hi-Air mixtures 


<l> 

Vcj 

(m/s) 

Tcj(K) 

Pcj(MPa) 

1.0 

1964.8 

2942 

1.58 


Table 5. Single-step reaction for stoichometric H 2 -Air mixtures 


<t) 

Reaction 

1.0 

H2 + ‘/2(02 + 3 . 76 N 2 ) -> H 2 O + 3.76/2N2 


An implicit scheme with a fixed time step of 1 nanosecond was first used during 
the initial comparison process. In order to validate the accuracy of the simulation run for 
each computational mesh size, the computed detonation velocity between each time step 
along the symmetry axis was computed using a first order finite difference method and 
compared against the C-J velocity. The C-J velocity is the velocity at which a steady state 
detonation wave propagates, which has been verified in shock tube experiments, as noted 
by Pickett and Davis [17]. 
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1. Computation of Detonation Velocity 

Detonation velocity was determined by noting the location of the shock front in 
the form of a pressure discontinuity along the symmetry axis of the cylinder 
configuration between each subsequent time step as illustrated in Figures 19 and 20. 



Figure 19. Pressure contour for Hi-Air mixture at 3 microsecond time step 



Figure 20. Pressure contour for H 2 -Air mixture at 4 microsecond time step 
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To obtain the wave speed, a five-point stencil first order finite difference method 
[18] as shown in Figure 21 is performed in one-dimension where x represents the spatial 
discretized location of the shock front and t represents the fixed time step pertaining to 
the shock front location. 


O 
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Figure 21. Schematic of the 5-point stencil method 


The detonation wave speed, is then obtained via: 


F, 


DET 


— x{t 2^ "1" Sx^ "^0 — 8x^t — 0 
12 At 


(16) 


2. Results from Initial Comparison 

Results of the normalized detonation velocity with C-J velocity based on cell 
mesh sizes of 0.25 mm (~ One ZND length), 0.5 mm (~ Two ZND lengths) and 1.0 mm 
(~ Four ZND lengths) for stoichiometric H 2 -Air mixtures against run distance along the 
symmetry axis of the PDE cylinder is plotted in Figure 22. Convergence plots for each of 
the simulation runs pertaining to each individual mesh size are included in Appendix B to 
verify that the resolved flow field solutions are bounded. 
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Figure 22. Plot of wave velocity/CJ velocity against run distance for different 

mesh sizes 


All 3 cell sizes indicate good agreement with the predicted C-J velocity towards 
the end of the PDE cylinder. The detonation waves decay towards the steady state C-J 
velocity at the end where the order of accuracy from the mean C-J velocity ranged from 
3% to 5%. Table 6 gives a statistical summary of the detonation wave velocities in terms 
of the standard deviation and the average normalized detonation wave speeds. This 
clearly implies that for simulation of a fully developed detonation wave propagating at C- 
J velocity, the cell sizes are clearly sufficient. However, the specified initial conditions 
were not enough to induce the ‘wave’ transient two-dimensional effects constituting the 
longitudinal and transverse wave interactions that have been observed in PDE detonation 
soot foil imprint experiments. Hence, further grid sensitivity studies were performed with 
the initial detonation wave propagating at temperature and pressure conditions specified 
by Zaev et al [8]. 
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Table 6. 


Set of evaluation parameters for different mesh sizes 


Mesh Size 
(mm) 

Average 
Detonation 
Wave Velocity 
(m/s) 

Standard Deviation 
of Detonation Wave 
Velocity (m/s) 

Average 

Detonation Wave 
Velocity/CJ 
Velocity 

0.25 (1 ZND Length) 

1892.58 

70.88 

0.96 

0.50 (2 ZND Lengths) 

1864.54 

77.87 

0.95 

1.00 (4 ZND Lengths) 

1898.81 

175.55 

0.97 


Also, the computational effort required to complete a simulation run for 0.25 mm 
took 8 days based on an implicit scheme with a fixed time step. To reduce and optimize 
the computational time, an implicit scheme with CFL number equal to one was adopted. 
In order to validate the implicit scheme based on CFL number, the same test case model 
using a mesh size of 0.25 mm was simulated using these two employed schemes and 
shown in Figure 23. The results indicate good agreement with each other with the 
numerical accuracy not compromised by adopting a different scheme. 



Figure 23. Plot of wave velocity/CJ velocity against run distance for different 

implicit time schemes 
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B. FURTHER GRH) SENSITIVITY ANALYSIS 


Further analysis was performed on cell sizes ranging from 0.25 mm (~ One ZND 
length), 0.0625 mm (~ One-quarter ZND length) and 0.05 mm (~ One-fifth ZND length). 
A plot of the normalized wave speeds along the symmetry axis with the C-J velocity as a 
function of the run distance of the cylinder is shown in Figure 24. Initial conditions for a 
starting detonation wave are specified at a pressure of 3.95 MPa and a temperature of 
2942 K. 


Plot of Wave Velocity/CJ Velocity against Run Distance for Different 

Mesh Size 


1.60 

- Wave Velocity/CJ Velocity (Mesh Size 0.125 mm) 



Run Distance(m) 


Figure 24. Plot of wave velocity/CJ velocity against run distance for different 

mesh sizes for Courant number of 1 


At a computational mesh size of 0.0625 mm, initial wave transients could be 

observed for an under-driven wave propagating at a sub-CJ speed and increase after a 

small run distance to an over-driven detonation wave before decaying and oscillating 

about the steady state C-J velocity. This clearly shows the transition to detonation or the 

DDT process. Decreasing the cell size from 0.0625 mm to 0.05 mm shows no change in 

the initial transients as well as the decay behavior to the steady state C-J velocity, which 
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allows us to determine that 0.0625 mm would be a practieal eell size to capture the 
detonation physics, as well as apply to other fuel-oxidizer mixtures. 

Further insight could also be gained by first analyzing the pressure and 
temperature contours as well as the histories for the given mesh size of 0.0625 mm at the 
mid-point of the cylinder configuration as indicated in Figures 25, 26, 27 and 28 
respectively. 

Behind the shock front, the detonation wave structure comprising of transverse 
wave interactions along with the longitudinal propagating wave front are observed, which 
are clearly multi-dimensional effects for both pressure and temperature contours. The thin 
shock front layer is propagating at a pressure state of 2.4 MPa that is higher than the C-J 
pressure of 1.58 MPa. Just behind the shock front, the detonation wave structure 
approximately 45 mm long comprising of the triple points are close coupled, which is 
fueling the detonation wave process followed by an expansion zone. Similarly, the H 2 -Air 
mixture reacts at an almost constant C-J temperature condition of 3000 K throughout the 
detonation wave structure zone. 

From the pressure history data, the pressure state for lead shock front relaxes to 
1.0 MPa before increasing back up to 2.0 MPa. The initial decrease is due to the 
interaction of the transverse wave propagating from the side walls and the symmetry axis 
of the cylinder. However, when the triple points intersect, the waves reinforce each other 
and the superposition process increases the localized pressure points. This promotes 
conditions that sustain the propagation of the detonation wave. 
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Figure 25. H 2 - Air pressure contour plot at 65 mm 
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Figure 26. H 2 - Air pressure history plot at 65 mm 
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Figure 27. H 2 - Air temperature contour plot at 65 mm 
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Figure 28. H 2 - Air temperature history plot at 65 mm 
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The detonation wave is sustained throughout by analyzing the contour and 
histories pertaining to conditions of pressure and temperature at the cylinder end. This is 
illustrated in Figures 29, 30, 31 and 32. The same detonation wave structure closely 
coupled to the shock front is observed with no change in the reaction zone length with a 
larger expansion zone. The lead shock front is propagating at a slightly lower pressure 
state of 2.3 MPa. As before, the pressure decreases before increasing back up to 2.4 MPa 
due to the interaction of the wave transients. The reaction zone conditions are 
approximately constant at 3000 K. 

One interesting feature to take note is when the detonation wave structure is 
closely analyzed, as shown in Figure 33, the detonation cell width, which characterizes 
the triple point intersection of the mach stem, incident shock and reflected shock is 
approximately 8 mm, which was taken to be the largest cell size since the detonation cell 
width is irregular in nature. This coincides with experimental observations as predicted in 
Figure 15. 
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Figure 29. H 2 - Air pressure contour plot at 122 mm 
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Figure 30. H 2 -Air pressure history plot at 122 mm 
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Figure 31. H 2 - Air temperature contour plot at 122 mm 
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Figure 32. H 2 - Air temperature history plot at 122 mm 



Figure 33. Close-up view of detonation wave structure 
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This concludes the grid sensitivity analysis and establishes that the cell size of 
0.0625 mm (~ one-quarter induetion length) is suffieient to eapture the detonation wave 
physics in a PDE cylinder eonfiguration. The same computational cell size will be used in 
subsequent simulations with a modified eonfiguration for analysis of axi-symmetrie 
toroidal implosions. 

C. IMPLOSION ANALYSIS 

Using the modified PDE eonfiguration illustrated in Figure 18, a detonation wave 
was assumed to propagate at the inlet with stagnation pressure and temperature equal to 
3.95 MPa and 2942 K with the detonation produets obtained from NASA CEA. The 
reactant conditions were speeified to be at an initial pressure of 1 atm and temperature of 
298 K. This was similar to the eonditions simulated for the grid sensitivity analysis. The 
eell size used was 0.0625 mm as determined previously, whieh was required to eapture 
the eomplex flow field governing the shock physics. The same test ease model for 
stoichiometrie reaetion of H 2 -Air mixture assuming a ehemieal kinetie reaetion setup of 
18 reaetions and 9 species listed in Table 2 was used for the implosion analysis. 

In order to determine whether the propagating detonation wave can be sustained 
by the implosion proeess, the detonation velocity was evaluated along the mid-axis of the 
inlet (whieh would be the center of the implosion event) and along the symmetry axis of 
the shock tube eonfiguration as shown in Figure 34. The resultant detonation velocity was 
compared against the C-J velocity to see if the wave ean persist as a detonation wave. 
Therefore, the eorresponding analysis are deeoupled into two portions, the first pertaining 
to the mid-axis inlet (implosion) and the seeond pertaining to the symmetry axis. 
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Figure 34. Evaluation of C-J wave speed for implosion study 


1. Implosion Along the Mid-axis Inlet 

Results of the normalized detonation wave velocity to C-J velocity versus run 

distance along the mid-axis of the inlet were plotted in Figure 35. From the graph, the 

detonation wave velocity decreases initially up to a run distance of 15 mm (~ 0.4 of the 

total radius) below C-J velocity as it travels along the radial direction towards the 

symmetry axis. This can be attributed to rarefaction effects as the detonation wave starts 

to expand from the inlet entrance. This can explicitly illustrated from the respective 

snapshot pressure and temperature contours in Figures 36 and 37 for the imploding wave 

at a run distance of 10 mm. From Figure 36, there is an expansion zone whereby the 

starting detonation wave decreases in pressure from specified inlet stagnation conditions 

along the inlet axis all the way to the reaction zone attached to the lead shock front that is 

propagating at the wave velocity. This reaction zone is a region of compression where the 

pressures start to build up fueled by the chemical reaction of the fuel-oxidizer mixture. 

Notice that this zone is localized only at the tip of the spherical wave front due to the 

nature of the toroidal imploding geometry and is a typical example of a shock focusing 

technique. There is a large high temperature region approximately equal to C-J 

temperature throughout the spherical shock wave except for the region immediately to the 

left and right of the inlet entrance whereby no communication exists between the flow 

field within the shock and reactant conditions. The observed high temperature region is 

supported by the heat release due to the reaction of the fuel-oxidizer mixture due to the 
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interaction of the multi-dimensional wave structure as the shock front implodes towards 
the symmetry axis. The pressures, temperatures as well as the kinetic chemical reaction 
rates are clearly driving the spherical wave front. 



Figure 35. Plot of wave velocity/CJ velocity against tube radius along mid-axis of 

inlet 
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Figure 36. Pressure contour for imploding geometry at run distance of 10 mm 

along inlet axis 
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Figure 37. Temperature contour for imploding geometry at run distance of 10 

mm along inlet axis 
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As the wave propagates beyond 15 mm towards the symmetry axis, the wave 
velocity starts to increase back upwards towards C-J velocity by virtue of the imploding 
geometry. This illustrates the fact that the detonation wave is continuously supported by 
the shock focusing process that can be determined from the C-J velocity. The pressure 
and temperature contours corresponding to a run distance of 20 mm are illustrated in 
Figures 38 and 39 respectively. The spherical wave front from the previous run distance 
of 10 mm has now transformed into little spherical wavelets coalescing into the leading 
edge of the detonation wave. This is due to little compression waves that are accelerated 
due to generation of heat release from the reaction of the fuel-oxidizer mixture. This heat 
release form high temperature regions that increases the sonic velocity of the gas 
products tries to catch up with the leading edge of the shock front. Also, extremely high 
localized pressure regions denoted by tiny red regions are also observed at the leading 
edge that is at least a factor of two greater than the C-J pressure. A look at the 
temperature contour reveals an interesting phenomenon, where transverse wave 
interactions are observed from these spherical wavelets at the protruding edge of the lead 
shock front extending for a distance of 10 mm along the mid-axis inlet as shown in the 
close-up view for the same temperature contour in Figure 40. There is also an extremely 
thin region separating the ambient temperature from the high temperature region - 
characteristic of a shock leading front. Further away from the mid-axis of the inlet on 
both sides, the thickness of this region increases, which indicates that the reactions are 
falling behind the shock front, and the detonation process is decaying due to expansion 
effects. 
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Figure 38. Pressure contour for imploding geometry at run distance of 20 mm 

along inlet axis 
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Figure 39. Temperature contour for imploding geometry at run distance of 20 

mm along inlet axis 
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Figure 40. Close-up view of temperature contour 

The pressure and temperature history plots along the mid-axis inlet are used to 
quantify the high pressures and temperatures driving the imploding wave as shown in 
Figures 41 and 42. The lead shoek front is propagating at a pressure of 2.5 MPa, whieh 
essentially represents an overdriven detonation wave. Immediately behind the shoek 
front, the pressure drops drastieally to 1.2 MPa followed by a constant pressure region 
approximately 8 mm thick. As for the temperature at the leading edge, it increases slowly 
from 2500 K in an oscillating manner to almost 2900 K, which equals approximately to 
CJ temperature before spiraling downwards to 2000 K. This spiraling behavior is due to a 
huge expansion effect where the pressure decreases from the constant pressure region to 
ambient pressure conditions. 
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As the wave continues to implode towards the symmetry axis, the localized high 
pressure region attached to the shock front expand radially outwards with the wave at a 
run distance of 38 mm, which is almost the radius of the PDE tube configuration as 
shown in Figure 43. The pressure history in Figure 44 verifies the fact that the pressure 
continuously increases beyond the C-J pressure to 3.5 MPa and this further reinforces and 
builds up the detonation wave speed. Similarly, the expansion wave region thickness 
increases as the high pressure relieves itself to 0.9 MPa before a sharp decrease to 
ambient pressure conditions. This indicates that a second shock front is formed in front 
of a set of supersonic jets of products leaving the inlet region and entering the detonation 
tube. This phenomenon will be treated robustly in future endeavors. 
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Figure 43. Pressure contour for imploding geometry at run distance of 38 mm 

along inlet axis 
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Figure 44. Pressure history of imploding wave at run distance of 38 mm along 

inlet axis 

Similar to the observed pressures, there is a radial expansion of the high 
temperature region corresponding to approximately C-J temperature outwards from the 
inlet mid-axis as illustrated in the pressure contour of Figure 45. The temperature history 
plot shows an almost constant temperature region of 2900 K extending from the shock 
front to a thickness of 20 mm, which is almost half the radius of the PDE tube. 

Analyzing the entire propagation of the toroidal wave front motion proves that 
even though the detonation wave initially loses its speed due to sudden expansion effects, 
the implosion process coupled with the reaction of the fuel-oxidizer mixture strengthens 
and increases the wave velocity as it propagates towards the symmetry axis. The pressure 
state just before the symmetry also shows that the initial wave builds up to an overdriven 
detonation shock wave propagating by a factor of 2.2 of the C-J pressure due to the 
implosion process. 
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Figure 45. Temperature contour for imploding geometry at run distance of 38 

mm along inlet axis 



Figure 46. Temperature history for imploding geometry at run distance of 38 

mm along inlet axis 
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2. Implosion Analysis Along Symmetry Axis 

Looking at the passage of the shock front along the symmetry axis, it starts out at 
a highly overdriven state with a wave velocity propagating by a factor of 1.8 times the C- 
J velocity. This is a highly transient state that will eventually decay to the steady state C-J 
velocity as the wave propagates along the symmetry axis towards the end of the PDE 
configuration. This is shown in Figure 47. Note that the shock wave passage starts out at 
a run distance of 23 mm along the symmetry axis due to the implosion process at the 
mid-axis inlet. The first two run distance points are also discarded due to the 
implementation of the five point stencil finite difference method used to determine the 
wave velocity. 



axis 


The snap shot pressure and temperature contours in Figures 48 and 49 denoting 
the point just after impact of the implosion process at a run distance of 23 mm along the 
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symmetry axis indicates extremely high localized pressure and temperature regions just 
aft of the shock front. The multi-dimensional effects of the wave interactions are 
explicitly observed in the temperature contour. The pressure and temperature history 
evaluated along the symmetry axis in Figures 50 and 51 allows for quantification of the 
pressure and temperature states arising from the collision of the two implosion shock 
waves. The localized pressure state along the symmetry axis increased to as high as 6.4 
MPa with the corresponding temperature at 3700 K. In comparison to the characteristic 
CJ pressure, this is on the order of at least 4 times greater and the detonation wave is 
highly sustained by the implosion. 
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Figure 48. Pressure contour for imploding geometry at run distance of 23 mm 

along symmetry axis 
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Figure 49. Temperature contour for imploding geometry at run distance of 23 

mm along symmetry axis 
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As the detonation wave front propagates along the symmetry axis to a run 
distance of 50 mm as shown in the pressure and temperature contours of Figures 52 and 
53, it catches up with the spherical wave front extending from the symmetry length to the 
wall of the tube configuration. The thin region behind the spherical wave is still clearly 
being driven by the reaction process as it stays attached to the shock front above C-J 
pressure and temperature while the region along the symmetry axis are propagating at 
extremely high overdriven pressures and temperatures fueled by the collision of the 
imploding shocks. The pressure and temperature history plots along the symmetry axis in 
Figures 54 and 55 verifies this with the lead shock front at a pressure and temperature 
state as high as 6.0 MPa and 4000 K respectively. 
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Figure 52. Pressure contour for imploding geometry at run distance of 50 mm 

along symmetry axis 
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Figure 53. Temperature contour for imploding geometry at run distance of 50 

mm along symmetry axis 
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Figure 55. Temperature history of imploding wave at run distance of 50 mm 

along symmetry axis 






At a run distance of 98 mm towards the end of the tube, the pressure and 
temperature contours of Figures 56 and 57 depiets the merging of the two separate waves 
into an almost planar wave front propagating at C-J veloeity. The loealized high pressure 
region stays highly attached to the shoek front followed by an expansion region. A triple 
point diseontinuity can be observed further upstream of the expansion region denoting the 
interaetion of the longitudinal and transverse wave from the eonfmed walls of the shoek 
tube. The pressure and temperature states along the symmetry axis pertaining to the 
leading edge of the shock front are still at 2.0 MPa and 3000 K respeetively from Figures 
58 and 59 above C-J conditions, which highlight the fact that a detonation wave still 
persists and will eontinue to propagate along the length of the tube. 
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Figure 56. Pressure contour for imploding geometry at run distance of 98 mm 

along symmetry axis 
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Figure 57. Temperature contour for imploding geometry at run distance of 98 

mm along symmetry axis 
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Figure 59. Temperature history of imploding wave at run distance of 98 mm 

along symmetry axis 


Analysis of the implosion process along both the mid-axis of the inlet and 
symmetry axis points to the fact that a detonation shock wave persists almost throughout 
the entire geometry and is sustained both by the imploding wave front. This highlights 
the potential of introducing a toroidal implosion process for PDE type applications and 
allows for larger diameters to be factored into such configuration designs such that 
shorter tube lengths can be factored in without the need to account for DDT run 
distances. Subsequently, a parametric study will be conducted to evaluate exactly the 
extent of the diameter increase such that a detonation shock wave can be sustained for 
such designs. 

D. PARAMETRIC STUDY 

Parametric studies were performed on tube diameters of 114.3 mm and 133.35 
mm respectively with the same inlet, reactant and boundary conditions specified for the 
toroidal implosion analysis. Note that the previous tube diameter of 76.2 mm from the 
implosion analysis was included as part of the parametric study. Current computational 
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resources could not support further exploration of tube diameters beyond 133.35 mm due 
to huge computational demands required and this could potentially take months to 
compute depending on the test case model. For the study, geometry of the shock tube as 
shown in Figure 18 remained the same with the exception of the diameter. The wave 
velocity was similarly evaluated along the mid-inlet and symmetry axis to determine if 
the detonation wave can be sustained. 

1. Parametric Analysis Along Mid-inlet Axis 

Figure 60 illustrates the detonation wave velocity as a function of the implosion 
shock tube radius for the different diameters specified along the mid-inlet axis. 



Figure 60. Plot of wave velocity/CJ Velocity against tube radius along inlet axis 
for different imploding shock tube diameters 


As the tube diameter increases, the detonation wave velocity starts out lower as 

under-driven detonations due to decreased confinement effects. An analysis of the 

pressure and temperature contours for diameters of 114.3 mm and 133.35 mm as the 

59 















imploding wave propagates mid-way through the radius of the tube shows lower pressure 
and temperature states behind the lead shoek front as eompared to the diameter of 76.2 
mm. Figures 61 and 62 illustrates the pressure and temperature eontours for 114.3 mm. 
Figures 63 and 63 illustrates the pressure and temperature eontours for 133.35 mm. 
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Figure 61. Pressure contour for imploding geometry at run distance of 28.75 mm 
for tube diameter of 114.30 mm along inlet axis 
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Figure 62. Temperature contour for imploding geometry at run distance of 28.75 
mm for tube diameter of 114.30 mm along inlet axis 
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Figure 63. Pressure contour for imploding geometry at run distance of 33.34 mm 
for tube diameter of 133.35 mm along inlet axis 
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Figure 64 . Temperature contour for imploding geometry at run distance of 33.34 
mm for tube diameter of 133.35 mm along inlet axis 


For the tube diameter pertaining to 114.30 mm, the phenomena corresponding to 
the coalescence of tiny spherical wavelets representing localized high pressure regions at 
the leading edge of the toroidal wave front could be observed similar to the implosion 
analysis. However the pressure states in these localized regions were lower in comparison 
due to geometrical differences in diameter. Similarly in comparison to the temperature 
contour for the implosion analysis, there is a thin separation region between the reactant 
conditions and the high temperature region at the toroidal front with an increasing 
thickness further away on both sides of the mid-axis inlet indicating a decaying 
detonation process due to expansion effects. These effects are however greater as 
compared to the test case model for 76.2 mm used for the implosion analysis and this 
again is a function of the diameter increase. This would account for the lower detonation 
velocities observed as the tube diameter is increased from 76.2 mm to 114.30 mm. 

For the case study of 133.35 mm, the pressure contour indicates a weak decaying 
wave propagating mid-way towards the symmetry axis. The pressure region behind the 
decaying shock front is almost at a uniformly low pressure with a magnitude of 
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approximately half of CJ pressure. However due to the geometry of the imploding tube, a 
seeond reflection wave from the head end of the tube can be observed interacting in a 
direction perpendicular to the weak wave propagation along the mid-inlet axis. This 
reflection wave represents a strong detonation wave with localized pressure regions just 
behind the toroidal wave front starting from the head end exceeding two times of C-J 
pressure. Likewise, as for the temperature contour, there exists a thick uniform region of 
separation between the reactant temperature and high temperature region along the 
toroidal wave front with an interacting reflection wave from the shock tube head end. The 
weaker pressure and temperature states behind the decaying shock front account for even 
lower detonation velocities as compared to the tube diameter corresponding to 114.30 
mm. 

As the wave implodes towards the symmetry axis, the process fueling the lead 
shock front is not supported and continues to weaken and decay due to expansion effects 
with the exception of the test case pertaining to 133.35 mm. The second reflected wave 
from the head end along the symmetry axis reaches interacts with the transverse wave 
propagation along the mid-inlet axis, which reinforces the lead shock front and results in 
a sharp increase in detonation velocity towards 0.8 of C-J velocity as shown in Figure 60. 
The respective pressure and temperature contours in Figures 65 to 66 and 67 to 68 
corresponding to the test case of 114.30 mm and 133.35 mm illustrate the observed wave 
propagation behavior as the shock front implodes towards the symmetry axis. 
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Figure 65. Pressure contour for imploding geometry at run distance of 52 mm 
for tube diameter of 114.30 mm along inlet axis 
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Figure 66. Temperature contour for imploding geometry at run distance of 52 
mm for tube diameter of 114.30 mm along inlet axis 
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Figure 67. Pressure contour for imploding geometry at run distance of 62 mm 
for tube diameter of 133.35 mm along inlet axis 
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Figure 68. Temperature contour for imploding geometry at run distance of 62 
mm for tube diameter of 133.35 mm along inlet axis 
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2. Parametric Analysis Along Symmetry Axis 

Along the symmetry axis, the detonation process is supported throughout as all 3 
test cases pertaining to different tube diameter geometries is observed to start from an 
overdriven detonation wave before decaying to a steady state wave propagating towards 
C-J velocity. A plot of the normalized wave velocity to C-J velocity against run distance 
along symmetry axis in Figure 69 depicts the observed wave propagation behavior. 


Plot of Wave Velocitv/CJ Velocity against Run Distance along Symmetry 

Axis 



Axial Distance (m) 


Figure 69. Plot of wave velocity/CJ Velocity against run distance along symmetry 
axis for different imploding shock tube diameters 


The test case pertaining to tube diameter of 114.3 mm starts out immediately at 60 
mm along the symmetry axis as the geometry of the wave front along the radial direction 
is such that at the point of implosion is almost planar throughout till the start point. The 
implosion process along the wave propagation drives the pressure and temperature along 
the symmetry plane to exceedingly high values to almost 50 times C-J pressure as shown 
in the pressure and temperature history of the starting wave at run distance of 60 mm as 
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illustrated in Figures 70 and 71. This eauses a sharp increase in C-J velocity to almost 1.7 
of C-J velocity, which causes the detonation wave to start out overdriven. There is also a 
second reflected wave from the head end of the tube, which is propagating behind the 
starting wave as shown in the pressure contour in Figure 72. 
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Figure 71. Temperature history of imploding wave at run distance of 60 mm 
along symmetry axis (diameter 114.3 mm) 
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Figure 72. Pressure contour for imploding geometry at run distance of 60 mm 

for tube diameter of 114.3 mm along symmetry axis 
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The detonation process is supported throughout the symmetry axis as both 
pressure and temperature states decay towards C-J conditions. Even as the wave approach 
towards the end of the tube at a run distance of 91 mm, the pressure and temperature 
states at the 1-D lead shock front along the symmetry axis as denoted by the respective 
temperature and pressure histories in Figures 73 and 74 are still above C-J conditions. 
This clearly indicates that the detonation propagation process is still supported. In 2-D, 
the pressure contour in Figure 75 at a run distance of 91 mm illustrates a non-planar 
shock front followed by a closely coupled high pressure region. In this region along the 
transverse direction, the pressure decreases further away from the symmetry axis before 
increasing again towards the wall of the shock tube. Although there is a decrease of 
pressure in the transverse direction, the states in this closely coupled are still above C-J 
conditions, which highlights the potential of increasing the diameter of the tube even 
further and still maintain sustenance of the detonation wave. 
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Figure 73. Pressure history of imploding wave at run distance of 91 mm along 

symmetry axis (diameter 114.3 mm) 
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Figure 74. Temperature history of imploding wave at run distance of 91 mm 
along symmetry axis (diameter 114.3 mm) 
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Figure 75. Pressure contour for imploding geometry at run distance of 91 mm 
for tube diameter of 114.3 mm along symmetry axis 
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For the test ease pertaining to the tube diameter of 133.35 mm, there are 2 
meehanisms sustaining the detonation wave propagation. The first mechanism is due to 
the reflection of the secondary wave from the head end of the shock tube while the 
second mechanism is due to the implosion process itself At implosion, the wave starts 
out overdriven with pressure and temperature states behind the lead shock front way 
above C-J conditions. For the specified geometry, the secondary reflected wave merges 
with the transverse imploded wave, which magnifies the pressures and temperatures 
behind the lead shock front. This can be observed in the pressure and temperature 
histories of the wave propagation evaluated along the symmetry axis at a run distance of 
30 mm in Figures 76 and 77, which reaches 32.5 MPa (~20 times C-J pressure) and 5000 
K (~ 1.7 times C-J temperature). The same 2-D pressure contour in Figure 78 illustrates 
an observed two distinct wave structure phenomena, a detonation wave that is 
propagating at super C-J velocity and a weak decaying wave that is remnant of the 
previous toroidal wave front that is slowing down due to expansion effects. 
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Figure 77. Temperature history of imploding wave at run distance of 30 mm 
along symmetry axis (diameter 133.35 mm) 
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Figure 78. Pressure contour for imploding geometry at run distance of 30 mm 
for tube diameter of 133.35 mm along symmetry axis 
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The super detonation wave travelling along the symmetry axis will eventually 
catch up and overtake the weak expansion toroidal wave front. This is illustrated in the 
pressure contour in Figure 79 as the wave approaches the end of the shock tube length at 
a run distance of 88 mm. The super detonation wave also expands towards the walls 
eventually replacing the weak expansion wave front as it propagates further along the 
symmetry axis. 
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Figure 79. Pressure contour for imploding geometry at run distance of 88 mm 
for tube diameter of 133.35 mm along symmetry axis 


The 1-D pressure and temperature histories of the same wave along the symmetry 
axis in Figures 80 and 81 still shows the detonation wave being supported even as it 
propagates towards the tube end. Here the C-J pressure and temperature states just behind 
the lead shock front are still at 4.0 MPa(~ 2.5 C-J pressure) and 3200 K (~ 1.1 C-J 
temperature), which will eventually decay to steady-state CJ velocity. 
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The parametric study highlights the immense potential of utilizing a shock 
focusing technique such as toroidal implosion to increase specific thrust in pulse 
detonation engines with no requirements to factor in tube length for DDT transitions as 
the diameter of the tube is further increased. Even when the diameter was increased to 
beyond 1.5 times the original length, the detonation process was still supported 
throughout as shown by the simulation results. The mechanisms supporting the 
detonation process was due to the implosion process itself as well as the secondary 
reflection wave from the tube head end. It would have been ideal to obtain the optimal 
tube diameter at which the implosion process fails to support the detonation process; 
however it was not possible to go beyond the tube diameter of 133.35 mm due to the 
limitation of the available computational resources. A two-fold increase in tube diameter 
results in a four-fold increase in required computational resource due to the number of 
mesh size. It is evident from the outcome of the parametric study that diameters way 
beyond 133.35 mm is possible to support the detonation process due to the simultaneous 
interaction and reinforcement of not only the implosion process itself, but also 
supplemented by the secondary reflection wave from the head end of the tube. 
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V. WAY FORWARD AND RECOMMENDATIONS 


Steps to conduct further parametric studies to determine the optimal diameter 
such that detonation process can be supported will be addressed in future endeavors. The 
specialized CFD solver from Metacomp Technologies requires extensive parallel node- 
to-node communications based on the computational cell size required to capture the 
complex physics governing the shock wave propagation process. This results in latency 
issues that contribute to substantial computational downtime whereby current 
supercomputing architecture on campus using infmiband interconnects is unable to 
support. Subsequent studies will be ported over to a high performance computing 
architecture off-campus with minimal latency and as such simulation time can be 
drastically reduced in order to allow for reasonable time for analysis of results. 

Additionally, it is desired not only to maximize the tube diameter based on the given 
geometry, but also the location of the inlet for the starting toroidal wave front. The location 
of the inlet from the head tube end will not only optimize the imploding detonation wave 
towards the symmetry axis but also optimizes the secondary reflected wave as it catches up 
with the toroidal wave front and aids as a supplementary mechanism to support the 
detonation wave propagation throughout the shock tube along the symmetry axis. 

Once the optimal configuration geometry has been established for imploding 
detonations for hydrogen-air, future endeavors will also address test cases for 
hydrocarbon fuels such as ethylene-air and JP 10-air. Moreover, optimal tube geometries 
can be further determined by using the practical computational cell size established for 
this thesis study as well as a set of viable reduced set of chemical kinetic reaction models 
specific to these fuel-oxidizer mixtures. Note that JPIO is a military qualified jet aviation 
fuel that due to its high energy density represents the best realizable fuel of choice for 
PDE applications. 

The future should focus on building an experimental test rig based on the 
imploding geometry to not only validate the simulation results, but in the process allow 
for a paradigm shift in future PDE designs that need not factor in multiple tube 
configurations and as such maximizing specific thrust. 
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VI. CONCLUSION 


A novel method to inerease specifie thrust in PDEs via imploding detonations has 
been investigated in this study. High fidelity simulations were required in order to capture 
the physics governing such complex flow fields. In this study, a practical computational 
mesh size on the order of one-quarter induction length was determined and deemed 
sufficient based on a reduced set of 18 reactions and 9 species chemical reaction model 
for hydrogen-air mixture. It was also found that the corresponding computed detonation 
cell width coincided with actual experimental data for the same fuel-oxidizer mixture. 

A set of simulations were carried out on an enhanced PDE design configuration 
utilizing implosion concepts and results indicate that the detonation wave was supported 
with pressure and temperature states persisting above CJ conditions behind the lead shock 
front throughout the entire geometry. Further parametric studies were conducted to 
evaluate and determine the extent of diameter increase for such design configurations and 
it was found that the detonation propagation process was still well supported even up to 
1.75 times the original diameter. This represents a three-fold increase in thrust per PDE 
tube, which roughly translates into a two-fold increase in specific thrust for similar total 
thrust PDE systems. However, there is a need to evaluate much bigger diameters. Due to 
the limitations of current parallel computing architecture, future endeavors will be 
addressed by using the Department of Defense (DoD) High Performance Computing 
systems that consists of a parallel computing platform with minimal latency in order 
determine the maximum diameter at which the detonation structure fails to persist in the 
specified geometry. 
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APPENDIX A: NASA CHEMICAL EQUILIBRIUM ANALYSIS 

OUTPUT 


NASA-GLENN CHEMICAL EQUILIBRIUM PROGRAM CEA2, FEBRUARY 5, 2004 

BY BONNIE MCBRIDE AND SANFORD GORDON 
REFS: NASA RP-1311, PART I, 1994 AND NASA RP-1311, PART II, 1996 

■k-k-k-k-k-k-k-k 


output massf 

prob case=00002142 det 

phi=l 

p (atm) =1 
output short 
output trace==le-10 


reac 

oxid Air t,k= 298 wt%= 100. 
fuel H2 t^k= 298 wt%= 100. 
end 

DETONATION PROPERTIES OF AN IDEAL REACTING GAS 

CASE = 00002142 

REACTANT 

TEMP 

K 

OXIDANT Air 

298.000 

FUEL H2 

298.000 

0/F= 34.29623 %FUEL= 2.833164 R,EQ.RATIO= 1.000000 PHI,EQ.RATIO= 

1.000000 


UNBURNED GAS 


PI, 

BAR 

1.0132 

Tl, 

K 

298.00 

HI, 

KJ/KG 

00 

1 

Ml, 

(1/n) 

21.008 

GAMMA1 

1.4015 

SON 

VELl,M/SEC 

406.6 


WT FRACTION 
(SEE NOTE) 

1.0000000 
1.0000000 


ENERGY 

KJ/KG-MOL 

-129.895 

-4.325 


BURNED GAS 
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P, BAR 

15.801 

T, K 

2943.63 

RHO, KG/CU M 

1.5500 0 

H, KJ/KG 

1333.33 

U, KJ/KG 

313.90 

G, KJ/KG 

-29833.7 

S, KJ/(KG)(K) 

10.5880 

M, (1/n) 

24.008 

(dLV/dLP)t 

-1.00954 

(dLV/dLT)p 

1.2062 

Cp, KJ/(KG)(K) 

3.3542 

GAMMAS 

1.1637 

SON VEL,M/SEC 

1089.2 

DETONATION PARAMETERS 

p/pi 

15.595 

T/Tl 

9.878 

M/Ml 

1.1428 

RHO/RHO1 

1.8042 

DET MACH NUMBER 

4.8335 

DET VEL,M/SEC 

1965.1 

MASS FRACTIONS 

*Ar 

1.2550-2 

*CO 

1.3093-4 

*C02 

2.6525-4 

COOH 

2.6149-9 

*H 

2.4760-4 

HCO 

1.0225-9 

HNCO 

2.395-10 

HNO 

3.1235-6 

HN02 

8.9839-7 

HNO 3 

1.028-10 

H02 

1.5250-5 

*H2 

2.6254-3 

HCOOH 

3.735-10 

H20 

2.2035-1 

H202 

2.3080-6 

*N 

1.1146-6 

*NH 

4.5040-7 

NH2 

2.7000-7 

NH3 

3.4969-7 

NH20H 

4.217-10 

*NO 

9.4235-3 

N02 

6.6339-6 

*N2 

7.2938-1 

N20 

2.7711-6 

*0 

1.3571-3 

*0H 

1.3522-2 

*02 

1.0117-2 

03 

6.2078-9 
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* THERMODYNAMIC PROPERTIES FITTED TO 20000.K 


NOTE. WEIGHT FRACTION OF FUEL IN TOTAL FUELS AND OF OXIDANT IN TOTAL 
OXIDANTS 
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APPENDIX B: CONVERGENCE PLOTS FOR RESIDUAL ERRORS 


A. GLOBAL RESIDUAL ERROR PLOT FOR MESH SIZE OF 0.25 MM 



B. INNER RESIDUAL ERROR PLOT FOR MESH SIZE OF 0.25 MM 
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GLOBAL RESIDUAL ERROR PLOT FOR MESH SIZE OF 0.50 MM 
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D. INNER RESIDUAL ERROR PLOT FOR MESH SIZE OF 0.50 MM 
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E, 


GLOBAL RESIDUAL ERROR PLOT FOR MESH SIZE OF 1.00 MM 
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G, 


GLOBAL RESIDUAL ERROR PLOT FOR MESH SIZE OF 0.0625 MM 
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I, 


GLOBAL RESIDUAL ERROR PLOT FOR MESH SIZE OF 0.05 MM 



J. INNER RESIDUAL ERROR PLOT FOR MESH SIZE OF 0.05 MM 
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K, 


GLOBAL RESIDUAL ERROR PLOT FOR IMPLOSION (76.2 MM) 



L. INNER RESIDUAL ERROR PLOT FOR IMPLOSION (76.2 MM) 
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M. GLOBAL RESIDUAL ERROR PLOT FOR IMPLOSION (114.30 MM) 



N. INNER RESIDUAL ERROR PLOT FOR IMPLOSION (114.30 MM) 
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GLOBAL RESIDUAL ERROR PLOT FOR IMPLOSION (133.35 MM) 



P. INNER RESIDUAL ERROR PLOT FOR IMPLOSION (133.35 MM) 
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